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I. INTRODUCTION 

The Theological properties of a suspension will depend, 
when the particles are non-spherical, on the orientation 
taken by the particles in response to the external flow. 
For a few particle shapes (e.g. the case of the ellipsoid 
0), equations for the rotation dynamics exist in closed 
form, and it is possible to determine the orientation dis- 
tribution of the particles in suspension. However, unless 
a mechanism for the achievement of a statistical equi- 
librium is introduced, the orientation distribution will 
depend on the state in which the suspension is prepared 
initially. In the case of microscopic particles, one such 
mechanism is provided by Brownian rotations It is 
still unclear whether inertia and interaction with other 
particles may contribute to the equilibration mechanism. 

An equilibrium distribution could be achieved alterna- 
tively by the presence of chaos in the rotation dynamics; 
unfortunately, the importance of chaos turns out to be 
small in most situations. In the case of a simple shear 
and axisymmetric particles, the particle motion is peri- 
odic Ij. This motion becomes aperiodic in the case of 
a time-dependent flow, but remains non-chaotic for ax- 
isymmetric particles Sj. Chaos arises in the motion of a 
triaxial ellipsoids in a simple shear but, depending on 
the axis ratios, large domains of initial conditions remain 
associated with regular orbits and to the absence of a 
uniquely deflned equilibrium distribution. Furthermore, 
for weak Brownian motion, the regular regions will act 
as attractors for the chaotic orbits and will provide the 
bulk of the orientation distribution. 

The equilibrium orientation distribution of a Brownian 
particle has been determined in various important limit 
regimes, depending on the value of the Peclet number Pe, 
defined as the ratio of the velocity gradient and the angu- 
lar diffusivity (which has the dimension of a frequency). 
The case of strong Brownian rotation was considered by 
Burgers |^ leading to an orientation distribution that 
in first approximation can be considered isotropic. A 
systematic perturbation theory in powers of Pe was in- 
troduced in ^6], allowing the calculation of the effective 
viscosity of dilute suspensions, for values of Pe up to 
20-30 0. 



More interesting is the case of weak Brownian motion, 
in which the form of the equilibrium distribution is de- 
termined by the structure of the orbits in orientation 
space, which in turn depends on the imposed flow. A 
technique for the determination of the equilibrium distri- 
bution of weakly Brownian particles, based on singular 
perturbation analysis of the diffusion equation in orien- 
tation space, was derived in |^ for the case of axisysm- 
metric particles in a simple shear. 

In the present paper, an alternative approach will be 
presented, based on the perturbative determination of 
the orbits in orientation space. This approach will appear 
to be appropriate in the case the flow is time-dependent, 
and, more in general, when analytical expressions for the 
unperturbed orbits are not available. For the sake of def- 
initeness, the dynamics of a small disk in the field of an 
oscillating strain flow will be analyzed, and its contribu- 
tion to the medium effective viscosity determined. This 
kind of flow could be obtained by means of a four-roll 
mill, as described in ^^J; more interestingly, as it will be 
illustrated in the next section, an oscillating pure strain 
is what is seen by a particle in the velocity field of a grav- 
ity wave. This turns out to have application to models of 
wave propagation in polar seas. In fact, under cold and 
windy conditions, high concentrations of millimeter size 
ice crystals, called frazil ice, are generated in the water, 
modifying the medium viscosity and leading to increased 
wave damping 9. lOj. (Given the crystal size, Pe is in 
this case typically very large). 

This paper is organized as follows. In the next section 
the orientation dynamics of a small disk will be analyzed 
in the absence of Brownian rotations. In Section III, the 
diffusion and drift across orbits in orientation space, pro- 
duced by a weak noise, will be calculated perturbatively. 
In Section IV the results will be applied to the calculation 
of the effective viscosity of a dilute small disk suspension. 
Section V is devoted to conclusions. 



II. UNPERTURBED ORIENTATION 
DYNAMICS 

Consider the motion of a particle in a velocity field 
U = {Ui, 1/2,0), which, at the particle position, has zero 



vorticity, and strain E = [VU + (VU)'^]/2 with compo- 
nents in the directions xi and X2'- 



E 



(1 + a) cosLot, 
(1 — a) sinLot, 



(1 — a) sintji 

-(1 + q) cos U!t 



(1) 



The interest in Eq. ^ is that it describes what is expe- 
rienced by a fluid element in a gravity wave. In fact, the 
velocity field of a small amplitude gravity wave, would 
read, choosing the X2 axis pointing downward from the 
unperturbed water surface at a;2 = llj: 



Ui = U[c-''''^ + e'=^(^^-2/0] sin(fcj;i - cut), 
U2 = U[e-''^^ - e'=(^2-2/i)] cos(fca;i - ujt), 



(2) 



where h is the water depth. From the potential nature of 
the flow, the vorticity of the flow O = [VU - (VU)T]/2 
is zero and the strain at the sea surface X2 — 0, has the 
form given by Eq. with e = kU, a = exp(— 2fc/i), and 
an initial phase different from zero for xi =/= 0. For small 
wave amplitudes, the displacement of a particle initially 
at the water surface will be small and the strain field 
experienced E(x(t),t) ~ E(x(0),t) will be given by Eq. 
(^. Expressions for the strain field in the form of Eq. 

can be shown to occur also for /i ^ 00, from the 
superposition of progressive and regressive waves, with 
e = fcf/+ > kU~ , a = U~/U~^, and the amplitudes 
of the two wave components. 

For a = 0, Eq. ^ describes a constant strain rotating 
with frequency lu/2 around the x^-axis, which, in the the 
gravity wave example, is associated with particle orbits 
that are perfectly circular JTIj. Transforming to the ro- 
tating reference frame leads to the new expression for the 
strain field 



E 



a sin 2ujt, 1 + a cos 2iLjt 
1 + a cos 2Lut, —a sin 2Lut 



(3) 



(the initial phase of the rotation has been chosen to pro- 
duce, at t = 0, a strain field with expanding direction at 
7r/4 with respect to the new xi axis; see Appendix A). In 
the rotating reference frame, an additional vorticity field 
is produced: 



n = ^ 1 

2 I -1 



(4) 



For a 7^ 0, this is a time-dependent planar flow, of the 
kind discussed in |^ , which is known to produce aperiodic 
behaviors in the particle orientation dynamics. 

The motion of a revolution ellipsoid, with symmetry 
axis identified by the versor p, in the presence of the 
strain and vorticity fields E and Jl, is described by the 
Jeffery's equations 1]: 



p = O p + G[E p- (p E.p)p]. 



(5) 



The parameter G gives the ellipsoid eccentricity, defined 
in terms of the particle aspect ratio r = a/h, where a 




FIG. 1: The coordinate system, 
rotating reference frame. 



The axes Xi are in the 



and h are respectively along and perpendicular to the 
symmetry axis, by means of the relation 



G = 



- 1 



1 



Introducing polar coordinates (see Fig. ^) and nor- 
malizing time and vorticity by the strain strength e: 
bj — > w/(— 2Ge) and t — > —Get, the Jeffery's equation 
(O, using Eqs. Q and Q, will take the form: 



?/. = -^ + /3(V',t), 



c=-i/3'(V,i)c 



(6) 



with c = i&nO, dot and prime indicating respectively 
d/At and d/d^, and: 



/3(?/;, t) = - cos2V' - acos{Aujt + 2V-'), 
I3'{ilj,t) = 2[sin2?/'4-asin(4cji-f 2'(/')] 



(7) 



(more details in Appendix A). 

Following 0, the orbits can be classified studying the 
Poincare map Pn{'4') = mod(-0(nT |'0),7r), with T = 
the period of (3, where ^{t\ip) obeys the first of Eq. © 
with ip{0\'ip) = ip. This eliminates the explicit time de- 
pendence from the dynamics, and will allow to isolate the 
slow, noise produced deviation between orbits, from the 
fast motion along them (see next sections). 

Some properties of the Poincare map can be obtained 
from inspection of Eqs. (|6I7|I . In particular, it is possible 
to see, from (3{tp, t) = /3{—ip, —t) and the form of Eq. lO, 
that the following relation holds: 



(8) 



and therefore the Poincare map is symmetric under the 
double reflection {ip, F„} {■jT — tp^n — P„} (see Fig. 
This has the consequence that fixed points, when present, 
would come in pairs located symmetrically around ip = 
tt/2 (see Fig. Hi). 

A fixed point in the Poincare map will be associated 
with a periodic tp and correspond to coherent orientation 
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FIG. 2: Symmetry of the Poincare map for the dynamics of 
Eqs. 16171 . Values of the parameters: uj = 1.4, a ~ 0.37, 
n = 1; G < (oblate ellipsoid). From Eq. ||HJ, one has that 
P_„(7r — %l))=Ti — Pn{tp) and therefore the plot is symmetric 
under reflection across the diagonal line Pn — tt ~ tp. 



of the particles. This regime is clearly produced by the 
aligning effect of strain on the particle dynamics, which 
becomes dominant in the small lu range. In the present 
case, Eqs. (I6I7|I appear to lead at most to a pair of fixed 
points, of which the the stable one is located at ip > 7r/2 
(see Fig. The stable fixed point tends to -0 = 37r/4 
at w = 0, corresponding to alignment of the long ellipsoid 
axis with the strain expanding direction. 

A transition to a coherent orientation regime is pre- 
dicted in the case of a deep gravity wave ^| at the 
crossover frequency Wc = 1 (no superposition of regres- 
sive and progressive components). Now, at — 1, the 
small wave amplitude approximation, leading to Eqs. lO 
[Tj) ceases to be valid (back to dimensional units, one 
would have for the particle displacement Ax ~ U /lOc ~ 
k~^). Nonetheless, the linearized theory provides a qual- 
itatively correct picture, as a transition to a coherent ori- 
entation regime is observed experimentally in wave tanks, 
although with a larger transition frequency ujc — 1-43 0. 
The presence of a coherent orientation regime appears to 
be preserved for a ^ with a crossover frequency Uc 
slowly decreasing as a ^ 1 (at a ~ 0.82, one has still 
LVc — 0.7). The decrease in the crossover frequency can 
be explained in terms of the destabilization of the fixed 
particle orientation in the rotating frame, by the oscillat- 
ing strain component of Eq. . 

The alternative regime of random particle orientation, 
is associated with |i/)(nT)| increasing monotonously with 
n, with Pni'4') generally aperiodic. In this case, from con- 
tinuity of 7/'(t|V'), Pn{'4>) will be topologically equivalent 
to an irrational rotation, and the sequence Pni'tp) origi- 
nating from a single tp will fill densely the interval [0, tt]. 



An ergodic property is then satisfied, i.e. it is possible to 
calculate averages over ijj as time averages. Furthermore, 
from Poincare recurrence, the Pni'ip) sequence will come 
arbitrarily close to the initial condition for some n. 

An important property is the following: if the or- 
bit starting from a certain tjj is approximately closed 
at t = nT, as shown in Fig. ^p, the same will occur 
with the orbits starting from any other initial condi- 
tion. In fact, if Pn{4') — "0 is small, the same will be 
true also for PniPm{tp)) — Pmi^) with m arbitrary. This 
is consequence of the dynamics of ijj being not chaotic 
(i.e. neighbouring trajectories do not separate asymptot- 
ically). Hence, exploiting the fact that, from ergodicity, 
Pm(V') fills densely the whole interval [0,7r], Pn{tp) — 
will be small for ip generic. 

Turning to the polar angle, if the orbit in ip is approxi- 
mately closed at nT, also c{nT\c, ■0) will come arbitrarily 
close to the initial condition c(0|c, ip) = c. Hence, to iden- 
tify approximately closed orbits, it is sufficient to look for 
recurrence of the Poincare map P„(0). To see why this 
property holds, the two of Eq. © can be integrated to 
give: 

, c{nT\c,^) 1^ dPni^) 

In = — in = — , 

c 2 dip 

and the condition c{nT\c,ip) ~ c will be satisfied pro- 
vided dPnjd'tp ~ 1. In the present case, the condition 
dPn/dip ~ 1 is satisfied provided Pni^p) — "0, i.e. if 
the orbit is approximately closed. The contrary would 
require dPn/dip to oscillate in -ip around dPn/difj = 1. 
However, if |Pn(0) — ipl < e for for some small e, the dif- 
ference dPn/dip — I could remain of 0(1) in intervals at 
most of length 0(e), in which one would have in conse- 
quence d'^Pn/d'ip'^ — 0{e~^). But this is prevented from 
smothness of the trajectories and of the functions f3 and 
/3'. 



III. THE EFFECT OF NOISE 

A. Coherent orientation regime 

Noise produces qualitatively different effects in the co- 
herent and in the random orientation regimes. In the 
coherent orientation regime, the main effect is smearing 
the transition to the random orientation regime. It is eas- 
ier to describe what happens at the transition for a — 0, 
where analytical expressions for ijjit) and c{t) are avail- 
able. When the crossover frequency = 1 is approached 
from above, i.e. from the random orientation regime, the 
rotation period for ip: T^, in the absence of noise, will 
tend to infinity like {uj^ — 1)^^ 0|. For lj < 1, the parti- 
cle is stuck at the stable fixed point = (cos^^ a; + 7r)/2 
and the period is by definition infinite. It turns out, 
that adding a small noise eliminates divergence of for 
UJ ujc, StS noise allows the particle to escape from the 
fixed point. 



a 







FIG. 3: Poincare map for the coherent orientation (a) and the random orientation regime (6) of an oblate ellipsoid in a shallow 
water wave with a ~ 0.37. Notice in case (a) the stable fixed point at ^ > it/2 and the unstable one at t/) < it/2. With prolate 
ellipsoids, the fixed points would have been exchanged. The value of n in the random orientation case (6) has been chosen to 
lead to approximately closed orbits. Notice that ip — tt/2 remains the best approximation to a fixed point (i.e. a closed orbit 
of period nT) . 



The time of escape from the fixed point can be ex- 
pressed in terms of the inverse of the probability for ■0 to 
reach the border of the basin of attraction for ^, which 
is tp = 7r/2; in other words: T'^ ~ < tt/2). The 
probabihty P(i/; < tt/2) could be roughly estimated, ap- 
proximating the dynamics of ip by the one of the Langevin 
equation obtained through linearization around Tp, in the 
presence of noise, of the first of Eq. (jSJ : 

dtp = 2{ip - i}) sin 2?/; At + dW. 

Here dW is the Brownian noise increment (Wiener pro- 
cess [13): {AW) = 0, {AW^) = At, and D is supposed 
small. This Langevin equation leads to the PDF (proba- 
bility density function) for if) [l3|: 

/ , X / 21 sin2V;| , , 
p[tp) — const, exp ^ — (ip — ^p) j . 

For small values of D, P{ip < tt/2) ~ p{tt/2) and there- 
fore: 

T.^exp(^i^(V^-./2)^). 

Thus, the rotation period is exponentially long in the 
inverse noise amplitude and the effect of Brownian rota- 
tions on the orientation distribution, which is governed 
by the stable fixed point of P„ (tp) , will vanish in the zero 
noise limit. 



B. Random orientation regime 

In the random orientation regime, the role of noise 
in determining an equilibrium orientation distribution is 



fundamental. If Brownian rotations were strictly zero, 
the evolution of the PDF p{ip, c; t) would be given by 
propagation along the unperturbed trajectories described 
by Eq. JBJ , which, from now on, will be identified by sub- 
script zero: 

piMtH), co{t\c, ^);t)= p{{p, c; 0) J(V;, c), (9) 

where J{4',c) = | det[(dV'07 dco)/(d'0, dc)]|~-^ is the Ja- 
cobian of the transformation {'ip,c} — > {ip^,c^}. This 
PDF is itself recurrent at the recurrence times ti — UiT, 
i — 1,2,..., for which P^.{{p) ~ {p, and therefore also 
co(ii|c, {p) ~c (see the end of last section). Hence, mem- 
ory of any initial PDF p{c,ip\Q) would be preserved at 
arbitrary large ti. p{c,tp]ti) ~ p(c, ■0; 0) and no relax- 
ation to equilibrium would be possible. 

Adding noise allows to reach statistical equilibrium in 
a time of the order of the inverse of the noise ampli- 
tude. Restricting to the discrete times nT, to make 
the process stationary, the equilibrium PDF for -0 is 
obtained from ergodicity and is the unique stationary 
PDF /9£(0) = const. |w + /3(V',0)|^^ The statistics of c, 
can then be described in terms of the conditional PDF 
p(c\ip) — p£;(?A, c) //?£;('(/'), where pe{iP,c) is the equilib- 
rium joint PDF at the instants t = nT. Notice that, 
from ergodicity of ip, it is sufficient to prescribe the form 
of piclip) at a single position -0. In fact, to obtain p(c|0) 
at '0 7^ "01 to lowest order in the noise, it is sufficient to 
propagate Eq. Q to t = nT and exploit the fact that 
PnW = mod(0o(»^7'|0), tt) is dense in [0, tt]. 

The first step to obtain a kinetic equation for p{c\^p), 
is to calculate the noise produced trajectory separation 
{c(t|c, 0) — Co{t\c, 'ip), ipitl'ip) — '0o(i|V')} at the recurrence 
times t = ti., where, choosing appropriately i,;, the differ- 
ences co{ti\c, Ip) — c and P^. (-0) — can be made small at 
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pleasure. However, since the conditioning in p(c|')/') is at 
-0 and not at ipitlip), it is then necessary to correct the 
first step and calculate the deviation c{t\c, ^p) — co{t\c, ip) 
at the time ti (equal to ti only in the zero noise limit), 
at which mod(V'(ii|V'): "') ^^nd ip are strictly equal. This 
means considering, instead of a Poincare map synchro- 
nized with the period of /?, the one synchronized with 
the rotation period in ip, i.e. with the crossing of tp by 

V'. 

Accounting for the effect of Brownian rotations, the 
Jeffery's equations will read (see Appendix B): 



dV'= hc^ + /3(^,i)]dt + i?V2gi(c)dW^^, 

-\f3'{^,t)c + Df{c)]At + D^h^{c)dWc, ^ ' 



Ac- 



where D has the meaning of a diffusion constant (in the 
present dimcnsionlcss units, D^^ ~ Pe), dWfc, with k = 
ipyC, are the Brownian increments: 



(dW^fe) = 0, {dWkdWj) = 6kjdt, 



(11) 



and the functions /, g and h are given by (see again 
Appendix B): 



fic) 



9{c) 



i(l + c^)(i+c^), 
c 2 

2^2 



+ 1 and h{c) (1 + c^) 



(12) 



The unperturbed orbits, as already mentioned, indicated 
by {^o,Co}, obey Eq. ©: 



Co = — 2^0*^0' 



(13) 



with Pq = p){ipo, t) and similar definition for /3q. For small 
noise, the correction can be determined as an expansions 
in powers of £)^/^: ^p = ^/jq + V'i/2 + + ■■• and similarly 
for c, with the initial condition tpki^) = Cfe(O) — for 
A: > 0. The lowest order correction is obtained from 
linearization of Eq. (I10|l around {"00; cq}: 



dV'i/2 = /3^V^i/2di + D^/^g--dW^, 
dci/2 = [2/3000-01/2 - i/3oCi/2]dt - 



D^h^dWc. 



(14) 



From Eq. l(TT|l and from hnearity of Eq. ifTHl (0i/2) = 
(ci/2) = 0; but "01/2 a-nd Ci/2 are not zero. The covariance 
equations obtained from Eq. (|14|l are: 



'^(V'?/2)=2/3^(V.?/2)+i?5o, 



dt 



('01/2C1/; 
dt("^l/2) = 



!/ - ^/5o(V'l/2Cl/2) + 2(3oCQ{tpl/^) 

4/3oco(0i/2Ci/2) - /3o(ci/2> + -D/io, 



(15) 



and lead to a diffusion contribution to the deviation. To 
obtain the drift contributions, it is necessary to consider 
the next order in the expansion of Eq. Hl()|l . and the 
result is: 



2/3o(0? 



/2/ 



^(ci) =-i/3^(ci)+2/3oco(0i), 

+ 2/3o(V'l/2Cl/2) +/3oCo(V'V2> 



(16) 




¥+¥i/2+¥i 



FIG. 4: Orbit behavior in the proximity of the recurrent 
point X — c}; X — z unperturbed orbit; w — y noisy orbit. 
The deviation between orbits is identified by c. 



The lowest order contributions to diffusion and drift are 
therefore both 0{D), as they should. Some simplifica- 
tions of Eqs. (|13|l . (|15|l and (|16|l . taking care of the sin- 
gularities of Eq. (|12|l at c = 0, are still possible and are 
illustrated in Appendix C. 

Once the noisy trajectory {tp{t\'tp),c{t\c,'tp)} has been 
calculated up to the recurrence time U, it is necessary to 
follow it back to the time ti at which mod(0(ti 10), tt) = tp 
and calculate the difference 

C = c(£i|c, 0) - Co(tj|c,'0) ~ c(fj|c,0) - C. 

This operation is equivalent to the procedure, implicit 
in 8], of subtracting from the deviation {0(ii|0) — 

1pQ{titp),c{ti\c,'lp) - Co(ii|c,0)} ~ {01/2 + 01,Ci/2 + Cl}, 

the component along the unperturbed Jeffery's orbit, and 
keeping only the part associated with percolation be- 
tween the orbits. The necessary operations are illustrated 
in Fig. ^ and it is assumed that the orbits can be pa- 
rameterized locally with = 0o(i|'0) (this is possible if 
is chosen away from turning points). The first step is 
to calculate the difference in c between noisy and unper- 
turbed orbits, at the azimuthal angle -0 = ■0(ti|0) where 
1nodi^p{ti\^p),^^) ~ -0 + 0'i/2 + "01: corresponding to the 
points y and z in Fig. ^ To 0{D), the value of c at z is 

C+C^ C-f C.0(0i/2 -I- 0l) -I- ^C000i/2' (17) 

where and c^^ give the rise in c along the unperturbed 
trajectory: 



dco , d^co 
= _ and c^^ = ^. 



(18) 



with d/dip the derivative along the unperturbed orbit: 



d 

d0^ 



1 

00 



d_ 

dt 



■00 



d_ 

dip 



Co 



d_ 

dc 



(19) 
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which is calculated a,t ip = tp. Combining Eqs. (|18I19() 
with Eq. inji: 



_ /3'co 
'^'A - 2(aj-/3o)' 



(/3o'-/3o)eo 



2/3oCQ 



(20) 



where /3 = 9t/3 and similar for /?'. To obtain c, it is 
necessary to correct the difference in c between y and z, 
i.e. Ci/2 + ci — c, for the contribution from the deviation 
between unperturbed orbits, which, in the present case, 
is (ci/2 + ci - c)(-0i/2 + V'i)cv.c, where 

c^c = 9c^/i9c = Cq ^c^. (21) 

Working again to 0{D): 

C = Ci/2 + Ci - C - -01/22^0(01/2 - c) , (22) 

corresponding to a time along the trajectory: 

k ^U + [lo + 0)]- Vi/2 + 0{D). (23) 

Using the relation (^1/2^1/2) = ^Co(V''i) (see Appendix 
C), together with Eqs. lfT7|) and lf2^ . the following re- 
sult for the diffusion and drift across Jeffery's orbits is 
obtained, to 0{D): 



1/2/ + 4^^^1/2) ~ 2cv,(V'l/2Cl/2), 



{C) = {Ci) + {c^C^c - 5 )('/'?/ 



(24) 



and, combining with Eqs. H2Q(I and (|21(l . the noise in- 
duced deviation between Jeffery orbits is fully deter- 
mined. 



IV. DETERMINATION OF THE ORIENTATION 
DISTRIBUTION 



The quantities (c^) and (c) allow to determine the noise 
contribution to orbit deviation, at the corrected recur- 
rence times U, at which mod ('0(ii|'0), tt) = ip. Both 
quantities (c^) and (c) are obtained from integrals along 
the orbits and it is expected that an averaging process 
takes place, with {c^)/ti and (c)/i,; tending to finite lim- 
its as —^ 00. Integrating numerically Eqs. H13|l and 
()15I16|I [or, more simply, Eq. (C2)], with the initial con- 
dition {V;o(O|-0),co(O|c,-0 )} = { -0,c} = {0,0} and then 
substituting, with Eqs. (|20I21|) . into leads to the 

result in Fig. El Self-averaging of {(?)/ti takes place also 
for relatively large values of the tolerance e, which iden- 
tifies recurrence, and therefore the sequence U — ti{e) 
i — 1,2..., through the condition |P^('!/5) — ■(/'| < e. It 
is thus possible to introduce effective drift and diffusion 
coefficients a and D: 

aic,iP) = lim iri(c), Dic,^) = lim t^^c^), (25) 




FIG. 5: 



2000 



Determination of the normalized difTusivity 



{'^1/2) /Dt for different values of the tolerance e entering the 
recurrence condition |P„(-i/)) — < e with {tp,c} = {0,0}. 
Values of the parameters uj — 1.4, a ~ 0.37. Thin line e = 0.4; 
heavy line e = 0.1; diamonds e — 0.01 (the diamonds identify 
the actual position of the recurrence times). 



which will describe the dynamics of the Poincare map 



c{k) - c(0) = aU + D^'^[W{t,) - Wm. 



(26) 



This is a discrete Langevin equation, in which W{t) is 
again the Wiener process, with ([IF(i) — VF(0)]^) = t. In 
the small D limit, the recurrence times ti can be treated 
as continuous on the scale of the relaxation to equilib- 
rium. It is then possible to obtain a Fokker-Planck equa- 
tion for the evolution of the PDF for c(i), which, by con- 
struction, is nothing else than /3(c|'0). Now, from Eq. 
(ESI: 

p{c\^- k) - p{c\^; 0) = [1 + 0{D^/^)]Udtp{c\i^; t)\t=o, 

and, to lowest order in Z), it is possible to disregard the 
difference between ti and ti in p; this is equivalent to 
substitute c{ti) — c(0) into the left hand side of Eq. H26|l . 
Taking the continuous limit, leads to the Langevin equa- 
tion dc = a&t + D^/'^AW , which is associated with the 
Fokker-Planck equation |1^] : 



d,{ap)^\dl{Dp), 



(27) 



and the notation f, indicating a slow time scale, is a re- 
minder that Eq. H27I) is meaningful only at timescales 
long with respect to ti — i^-i. Slow variation of the flow 
parameters entering Eq. (|10|l would lead to dependence 
of the coefficients in Eq. H27|) on the slow time i. As 
in 1^], the fact that both a and D depend linearly in D 



7 




FIG. 6: Comparison of the PDF p{c\tp) calculated using for 
d and D different values of ti [see Eq. 1251 1. Values of the 
parameters: lj = 1.4, a = 0, tp = 0. Heavy line: ti = ST 
e = 0.4- stars: U = 20T, e = 0.1; thin line: Leal & Hinch 
theory Q]. 



implies that the equihbrium PDF is independent of the 
noise ampHtude. 

As statistical equilibrium will be achieved on the time 
scale D^^, in order for the approach to be meaningful, 
it is necessary that the ti used to define a and D satisfy 
Dti <C 1. Actually, excellent convergence is obtained 
already for ti rather small; in the case of Fig. at 
ti = 20T, corresponding to e = 0.1. Notice that the 
case considered in Fig. |H| which is identical to the deep 
water wave regime considered in [l^ , can be mapped to a 
constant simple shear by a redefinition of the eccentricity 
G. [In this case, the aperiodicity of -P^('0) originates not 
from the dynamics, but from the sampling time T and the 
rotation period Tr = (w^ — 1)~2 being incommensurate]. 
The PDF p{c\il>) can then be compared with the analyti- 
cal result from the theory of Leal and Hinch [the function 
/(C) in Eq.(17) of their paper] As can be seen from 
Fig. El the two approaches give indistinguishable results 
already for ti = 20T, e = 0.1. Similar convergence to the 
limit result is observed for a > 0, when comparison with 
the theory of Leal and Hinch is not possible. 

Knowledge of the PDF p{c\^p) allows determination of 
the effective viscosity of a dilute disk suspension in the 
oscillating strain field of Eq. . The viscous stress for 
a suspension of axisymmetric ellipsoids reads |^ Il4j , in- 
dicating with n and respectively, the solvent viscosity 
and the suspended phase volume fraction: 



(T = 2^E + 2^${2^(pppp) : E 
+2B[{pp) • E + E • (pp>] + CE}, 



(28) 



where, in the present time dependent situation, the av- 
erages are intended over orientation and time. The coef- 



ficients A — C depend on the particle geometry [l 

5 ^^'^ I B - ^ 64 1 
37rr Ott^ ' Snr Ott^ 2 

and C 



A 



8 ^ 128 
Sirr 97r^ ' 

with r the particle aspect ratio, supposed small. These 
expressions correct to 0(1), similar ones derived in j^. 
From the stress cr, the effective viscosity jl can be calcu- 
lated in terms of the viscous dissipation in the suspension: 



1 cr : E 



^{1 + K<^)fi, 



(29) 



2E : E 

where K is called the reduced viscosity. Expressing the 
versor p in function of the angles "0 and 0, and using Eqs. 
and 

K = A{s\n^ 9 sin^ 2^-) + 2B{sm^ 9) + C. (30) 

As in 0, the average over orientation is split into parts 
along and transverse to the orbit. In the present situa- 
tion, however, evaluation of the average along the orbit is 
slightly more delicate than in the time-independent case. 
At a generic time t the average of a function f{tp,c) will 
be: 



1 

(/)w = -E 

i—l 



dcp{c\^)f{i^{t + zT|^), c{t + iT\c, ^)), 



where, from ergodicity, memory of the initial condition 
-0 is lost for n oo. Carrying on the average over one 
period, which, in the case of a wave, from Eq. 0, is 
equivalent also to a space average, leads to the average 
along an orbit: 

TlT 



dcp(clV') / dtf{mi^),c{t\c,i;)). (31) 



Evaluating Eq. (|30|l with Eq. (|31|l leads to the values of 
the reduced viscosity shown in Fig. [7| The calculation 
has been carried on, using as recurrence point, tp = 0; 
the value of the particle aspect ratio has been chosen 
consistent with frazil ice measurements Q. The same 
qualitative regime observed for a = is reproduced here, 
namely, a dip in the reduced viscosity at the crossover 
from the coherent rotation regime to the random orien- 
tation one [i3] . 

In the case of gravity waves, the coherent regime would 
always be associated with high amplitude waves, corre- 
sponding to small values of the normalized frequency uj. 
The reduced viscosity has been calculated in this range 
from Eq. H30|l . fixing 9 — tt/2 and integrating the equa- 
tion for ip, in the absence of noise, with initial condition 
at the fixed point. 

V. CONCLUSIONS 



The numerical evaluation of the rheological properties 
of a suspension of particles that are weakly Brownian is 
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I ^ ^ ^ 1 

0.5 1 1.5 2 



FIG. 7: Reduced viscosity, averaged over a period, for a 
suspension of disk-like particles with aspect ratio r = 0.045. 
In the three cases: (a) a ~ 0.61; (b) a ~ 0.37; (c) a = 0. 



perimental data on the sea ice effective viscosity from 
wave tanks, in which the deep water condition is at the 
most only approximately satisfied, as a test case of what 
happens in the open sea |^ ^| . 

A natural extension of the techniques illustrated could 
be the treatment of higher numbers of degrees of freedom. 
An interesting example is the triaxial ellipsoid in a simple 
shear considered in In this case, the angle 9 would 
be replaced by the pair {9, 0} with </> the rotation around 
the axis p. An analysis in the whole phase domain would 
require, however, consideration of the transition region 
from the regular orbits, in which diffusion is dominated 
by Brownian rotation, to the chaos dominated stochastic 
region. 
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faced with difficulties associated with the long integra- 
tion times necessary to achieve statistical equilibrium. 
Analytical techniques for the calculation of the cumula- 
tive effect of the Brownian noise on the dynamics are 
therefore necessary. The technique presented in this pa- 
per can be seen as a multiple time scale analysis [T^ 
in which the stochastic dynamics is pushed to the slow 
scale, while the local strain and vorticity are treated as 
fast variables. For the periodic flows considered in this 
paper, the effective drift and diffusivity coefficients are 
obtained integrating the fast (and aperiodic) orientation 
dynamics up to the first recurrence time, at which the ap- 
proximation of a closed orbit is considered good enough. 
Slow variations would be accounted for, sampling the al- 
most closed trajectory segments in appropriate way along 
the particle orbit, and would lead to effective drift and 
diffusivity coefficients depending on the slow time. Once 
the effective drift and diffusivity were available, a Monte 
Carlo, for the determination of the rheological properties 
of a suspension, would be carried on at the slow time 
scale [this would be associated formally with integration 
of the Fokker-Planck equation H27|) ]. 

Application of these techniques to the dynamics of a 
thin disk suspension in gravity waves, suggest that qual- 
itative behaviors in deep water, accounted for theoreti- 
cally in 01 and observed experimentally in (Qj , should be 
preserved in the shallow water regime. In particular, a 
transition from a coherent rotation regime for large am- 
plitude waves to a random orientation one, marked by a 
deep minimum in the medium effective viscosity, should 
continue to be present. Away from this regime, in the 
random orientation regime, the numerical values of the 
effective viscosity appear to be only weakly dependent 
on the water depth. This strengthens confidence in ex- 



APPENDIX A: ORIENTATION DYNAMICS IN 
ROTATING REFERENCE FRAME 

The strain matrix in the laboratory frame, is provided 
by Eq. Passing to the rotating frame: 



R E 



where 



R 



l + acos2ti;t, — asin2cL;t 
— asin2a;t, — 1 — Q!Cos2cjt 



cos(jjt/2, smLot/2 
-sin(jjt/2, cosLot/2 



is the matrix for an angle —Lotjl rotation. In term of 
components, Ui and Vi = RijUj are the components of a 
vector in the laboratory and the rotating reference frame. 
In the rotating frame, the fluid is seen rotating like vf. 



2 



0, 1 
-1, 



and fl is the vorticity of the fluid measured in the ro- 
tating frame. One more rotation by 7r/4 produces Eq. 

m 



R 



1 + acos2u!t, 



-a sin 2ujt 



-asin2u;t, — 1 — acos2a;< 



a sin 2ujt, 1 + a cos 2Lut 
l + acos2ujt, —asm2Lut 



where 



9 



As the next step, pass to adimensional variables: 



t = -Get, 



2Ge' 



E = e^^E. 



This choice guarantees that, in the case of oblate ellip- 
soids, the signs of the normalized times and frequency 
are preserved. Substituting into the Jeffery's equation 
© gives: 



I )p-[E.p-(p.E.p)p]. (Al) 



Introducing polar coordinates p — (sin 9 cos ip, sin 9 sin -0, 
cos 61) and using Eq. ||2J), leads to the expressions: 

E.p = sin0(' + + (A2) 

\ cos ip + a cos[tp + Aiot) I ^ ' 



and 



p • E • p = sin^ 6'[sin 2i}) + a sva{2il) + 4a)i)]. (A3) 

The Jeffery's equation can now be written in components. 
Starting from 9, using Eqs. (A1-A3): 

P3, — ~ sin 99 — cos 9 sin^ 9 [sin 2iIj + a sm{2ip + iujt)] , 

(A4) 

which leads to the second of Eqs. Ht)l7|) . Passing to the 
equation for ip: 

P2 = sin 9 cos ipip + cos 9 sin ip9 — — (2> sin 9 cos ip 
— sin9[cos'ip + acos{i/j + Ailit)] + sin'^ 9[sm2ip 
+a sin(2-0 + 4tZ>t)] sin ip, 

from which, using Eq. (A4): 

cos ipip = —ijj cos Ip — cos 0(1 — 2 sin^ ip) 
-|-a[sin(2^ + Aujt) auiip — cos{ip + 4:Lut)], 

and, after little algebra: 

i/j = —Lu — COS 2ip — a cos{2ip + Aujt) , 

which is the first of Eqs. l|BI7|) . 



isotropic solution p{ip,9) = ^sin6']. One finds the in- 
crements for Ip and 9 produced by Brownian rotation in 
the time interval At [l^: 

dip =\ sin 9\-^dW^, A9 =]^cot9At + dW0, 

where dVFfc, k — ip,9 are the Brownian increments 

{dWk) = 0, (dWjdWk) = S.kdt. 

Changing then variables from to c and using Ito's 
lemma, one finds: 

dr 1 d^r 

= -(1 + + c')dt + (1 + C^)dWg, 

c 2 

and, using the expression sin6' = c(l + c^)^^/^ in dip, Eq. 
(|10|l is finally obtained. 



APPENDIX C: ALTERNATIVE FORM OF THE 
PERTURBED ORBIT EQUATION 



Equations (|15I16I) can be simplified, and the singular- 
ity in 9 = produced by the noise term in the first of 
Eq. (|1U|1 eliminated, by the change of variables: 

yi = C^(V'?/2>> y2 = Co(-01/2Cl/2>, 
2/3 = (Ci/2), 2/4=Co(ci), y5=Co(0i). 

In the new variables, Eqs. (|15ll()l) take the form: 



yi = P'oVi + Dg, 
2/2 = 2/3oyi, 
< 2/3 = 4/3oy2 - /3^,2/3 + Dh, 

iji = -Pi'oVi + 2poy2 + P'oVi + 2/3o2/5 + Df, 
12/5 = -2/3oyi, 



(CI) 



APPENDIX B: NOISE TERM DETERMINATION 

The noise term to add in Eq. © can be obtained di- 
rectly from the diffusion equation obeyed for zero flow 
by the orientation PDF in the variables {ip,c}. Alter- 
natively, one may consider the diffusion operator in the 
variables {ip, 9}: 

^ sin2 0ai/;2 + sm9d9''''' 89' 

and determine the stochastic process leading to the 
Fokker-Planck equation V'^p{ip,9) = [which has the 



where, from Eq. (fT^ . g — 1 + Cq, h = {\ + Cg)^ and 
/ = (l-|-Co)(5 -l-Co). Comparing the equations for 2/2 and 
2/5, one sees that, thanks to the initial condition yfe(O) = 
0, 2/2 = -2/5 and then (-01/221/2) = -co(0i); thus the 
equation for {ipi) can be eliminated from Equation 
(CI) can then be further simplified to: 

{2/1 = P'qVi + Dg, 
2/2 = 2/3o2/i, .^2) 
2/3 = 4/302/2 - P'm + Dh, 
2/4 = -/3^2/4 + /3o2/i + DJ. 
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